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ABSTRACT 


\ 

Finite  amplitude  internal  solitary  waves  in  a  stratified  fluid  are 
computed  numerically  as  solutions  to  a  version  of  Long's  equation  . 
Newton's  method  is  used  to  linearise  the  two  dimensional  nonlinear  elliptic 
equation  and  numerical  continuation  techniques,  both  using  the  wave  speed  and 
a  pseudo-arclength  parameter,  are  used  to  trace  out  solution  branches 
efficiently.  Numerical  results  for  the  'tanh'  density  profile  are  presented 
here  for  various  depths  of  the  fluid.  For  shallow  depths,  solutions  for  a 
fixed  wave  speed  are  not  unique  and  bore-like  solutions  with  large  amplitude 
have  been  found.  In  the  deep  water  case,  excellent  agreement  is  obtained  with 
the  experimental  data  of  Davis  and  Acrivos  4^1  whereas  traditional  weakly 
nonlinear  analysis  fails  to  produce  agreement  in  the  large  amplitude  regime. 


*This  work  was  supported  by  Dynamics  Technology  Inc.,  Torrance,  California 
and  by  ONR  Grant  N00014-80-0076  under  subcontract  from  Florida  State 
University. 
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1.  Introduction 

Large  amplitude  internal  waves  which  maintain  their  shape  over  long 
distance  of  propagation  (hence  the  name  ' solitary')  are  often  observed  in 
nature  and  have  recently  been  suggested  to  occur  in  the  Jovian  atmosphere  in 

the  form  of  the  Great  Red  Spots  [18].  These  waves  can  also  be  generated  in 

the  laboratory  rather  easily  [17] ,  which  suggests  that  they  can  be  excited 

under  many  circumstances  in  natural  systems.  Most  of  the  existing  theories 

(  [5],  [6],  [141#  [20]  )#  with  a  few  exception$(  [2]),  are  limited  to  the 
case  of  small  amplitude  and  long  wavelengths,  due  to  analytical  difficulties. 
The  small  amplitude  assumption  allows  the  solution  to  be  expressed  in  an 
asymptotic  expansion  in  terms  of  the  amplitude  and  the  nonlinear  governing 
equations  to  be  solved  by  perturbation  techniques.  The  long  wavelength 
assuatption  reduces  the  governing  two  dimensional  partial  differential  equation 
into  ordinary  differential  equations.  In  applications#  results  from  these 
weakly  nonlinear  theories  are  often  extrapolated  into  the  large  amplitude  and 
short  wavelength  domains.  While  it  is  remarkable  that  such  extrapolations 
agree  quite  favorably  with  experimental  data  even  for  moderately  large  wave 
amplitudes  and  short  wavelengths#  notable  discrepancies  exist  both 
quantitatively  and  qualitatively  in  many  situations. 

In  this  paper#  we  tackle  numerically  the  original  nonlinear  partial 
differential  equation  directly  without  making  either  the  small  amplitude 
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assumption  or  the  long  wavelength  assumption.  Large  amplitude  waves  with 
short  wavelengths  have  been  computed  together  with  the  small  amplitude  and 
long  wavelength  ones.  The  nonlinear  equation  is  linearized  by  Newton's  method 
and  the  resulting  linear  equations  are  discretized  by  standard  finite 
difference  methods.  In  the  deep  fluid  cases#  the  waves  only  occupy  a  small 
portion  of  the  computational  domain#  and  use  of  uniform  mesh  spacing  would  be 
wasteful.  Our  solution  is  to  use  an  analytical  transformation  of  the  domain 
so  that  a  uniform  mesh  in  the  transformed  domain  gives  adequate  resolution 
without  too  much  waste.  The  governing  equation  contains  a  physical  parameter 
the  wave  speed#  and  we  are  interested  in  computing  many  different  waves 
corresponding  to  different  wave  speeds.  To  make  the  computational  process 
more  efficient#  we  employ  continuation  techniques  (with  respect  to  the  wave 
speed)  whereby  known  solution  at  a  given  wave  speed  can  provide  extremely  good 
initial  guess  for  Newton's  method  at  a  nearby  wave  speed.  Moreover#  in  a 
shallow  fluid  case  we  found  nonunique  solutions  for  a  given  wave  speed.  In 
order  to  provide  a  systematic  way  to  trace  some  of  these  nonunique  solutions# 
an  arclength  continuation  technique  [12]  is  employed. 

In  Section  2#  the  governing  partial  differential  equations  are  discussed# 
and  in  Section  3  #  a  description  of  Newton's  method  and  the  continuation 

techniques  is  given.  Details  of  numerical  implementation  are  described  briefly 
in  Section  4#  and  numerical  results  and  comparison  with  experimental  data  are 
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presented  in  Section  5.  The  conclusion  is  in  Section  6. 

This  paper  is  a  companion  to  the  paper  by  Tong,  Chan  and  Kubota  [24]. 
There  we  focus  on  more  mathematical  issues  such  as  proving  existence  of  large 
amplitude  solutions  and  relating  the  dependence  of  wave  amplitude  on  wave 
speed  to  the  density  profile  of  the  fluid  medium.  For  a  more  thorough 
discussion  of  the  fluid  dynamical  aspects  of  the  waves  computed  in  the  current 
paper  and  other  background  information*  the  reader  is  referred  to  the 
companion  paper.  The  present  paper  focuses  on  the  numerical  techniques  used 
to  computed  the  large  amplitude  waves.  It  is  hoped  that  the  numerical 
techniques  are  of  more  general  interest  and  can  find  applications  in  other 
areas. 


2.  Governing  Equations 

Consider  a  stably  stratified  fluid,  the  ocean  for  example,  which  is 
perturbed  and  set  into  motion  (Figure  2-1),  The  interaction  of  the  density 
variations  and  the  force  of  gravity  can  give  rise  to  many  interesting  wave 
phenomena.  Internal  solitary  waves  are  a  well-known  example  [9,  10] .  These 
are  waves  that  propagate  inside  the  fluid,  rather  than  on  the  surface,  and  can 
travel  long  distances  without  changing  their  shape  or  losing  their  energy. 

The  equation  that  governs  the  steady  state  propagation  and  shape  of  these 
internal  waves  can  be  derived  from  the  equations  of  motion  of  an 
incompressible,  inviscid,  and  non-diffusive  fluid  of  variable 
density  [16,  24],  which,  in  terms  of  dimensionless  variables  for  the  case  of 
uniform  mean  flow,  is: 

An  +  XuF' (y+u)/ (l-oF(y+u) )  * 

0F'(y+u)(u*  +  uj  +  2  ny)/(l-oF(y+u))/2,  (2.1) 

with  boundary  conditions: 

u  *  0  on  the  surface  and  bottom  of  fluid, 
u  ->  0  as  I x I  ->  ®, 


where 
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Figure  2-1:  Internal  Solitary  Wave  Environment 


fluid  surface 


direction  of 


propogation 


p(y)  — 
density 


internal  wave 


"  u  is  the  perturbation  stream  function, 

-  X  is  a  nondimensional  quantity  inversely  proportional  to  the  square 
of  the  wave  speed, 

-  p  the  density  is  given  by: 

P  *  Pq(1  “  oF(y+u) ) ,  (2.2) 

where  F  is  a  given  nonlinear  function  that  specifies  the  density 
stratification  profile;  an  often  used  profile  is  the  r tanh'-prof ile : 

F(y)  =  tanh(y), 

-  a  can  be  interpreted  as  the  relative  density  change  across  the 
entire  stratified  layer. 

Equation  (2.1)  is  seldom  solved  in  this  complicated  and  fully  nonlinear 
form.  Often  the  relative  density  stratification  is  very  small  (e.g.  in  the 
ocean),  so  that  o  <<  1  and  the  so-called  Boussinesq  approximation  is  made. 
Equation  (2.1)  is  then  reduced  to  the  simpler  form,  to  which  we  shall  refer  to 
as  Long's  Equation  [16]: 

Au  +  XF'(y+u)  u  =  0.  (2.3) 

Equation  (2.3)  is  the  equation  we  will  be  dealing  with  in  this  paper. 
For  the  purpose  of  describing  the  numerical  methods  in  the  next  section,  it  is 


convenient  to  write  Equation  (2.3)  formally  as 


G(  u  ,  X  )  =  0  (2.4) 

to  show  the  dependence  of  the  perturbation  stream  function  u  on  the  parameter 
X.  Equation  (2.4)  is  in  the  form  of  a  nonlinear  elliptic  eigenvalue 
problem  [4f  12]. 


4 


\ 

1 


-  8  - 

3*  Newton's  Method  end  Continuation  Techniques 
3.1  Newton's  Method 

Solving  a  nonlinear  eigenvalue  problem  G(u,X)  =  0  can  take  on  at  least 
two  different  meanings: 

1.  Solve  for  u(X)  given  a  specific  value  for  X. 

2.  Determine  the  dependence  of  u  on  X* 

In  this  paper,  we  adopt  the  second  meaning.  Clearly#  if  we  have  a  good 
procedure  for  solving  Problem  (2),  we  can  also  use  it  to  solve  Problem  (1) 
effectively. 

The  usual  solution  procedure  for  nonlinear  problems  involves 
linearization  and  application  of  some  version  of  Newton's  method. 

Newton's  Method  : 

Given  a  value  of  X  and  an  initial  guess  u^  for  the  solution  u(X) #  we 
perform  the  following  steps  repeatedly  until  satisfied  : 


G*  5u*  =  “  G(u*#X) 
u 

(3.1) 

ni+1  =  a1  +  Su1. 

(3.2) 

In  the  above#  subscripts  denote  partial  derivatives  and  so  G^  denotes  the 
Jacobian  of  the  operator  G  (with  respect  to  u) .  This  procedure  will  generally 
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converge  quadratically  when  it  does  converge.  However,  as  is  well  known,  in 
many  instances  it  will  fail  to  converge  when  the  initial  gness  is  not  'close' 
to  the  true  solution. 

3.2  Natural  Continuation 

A  plausible  procedure  for  overcoming  this  convergence  difficulty  and  also 
for  determining  the  dependence  of  u  on  X  is  to  start  at  a  known  solution 
(u0,XQ)  on  the  solution  curve  and  use  it  as  initial  guess  for  a  Newton-type 
iteration  to  find  the  solution  for  a  neighboring  point  on  the  solution  curve 
with  X  close  to  X^  and  repeating  the  procedure.  We  can  improve  on  this 
further  by  computing  the  'slope'  u^  at  a  known  solution  and  use  it  to  get  a 
better  initial  guess  for  the  next  value  of  X  in  a  predictor-corrector  fashion. 
We  shall  call  this  the  Natural  Continuation  Procedure  because  it  corresponds 
to  parametrizing  the  solution  curve  by  X,  the  naturally  occurring  parameter. 

Natural  Continuation  Procedure : 

Given  a  known  solution  want  to  compute  the  solution  at  a 

nearby  value  of  X. 

1.  First  compute  the  'slope'  u^  at  (uq,Xq)  from 

Gu  ■  -  v  <3*3> 


2.  Perform  an  Euler  predictor  step: 
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u  *  o0  +  (X  -  XQ) . 

3.  Use  n°  as  initial  gness  in  Newton's  method  for  G(u.X)*0: 
G1  (Ui+1  -  u1)  -  -  G(ui,X) 

B 

until  convergence. 


(3.4) 


(3.5) 


4.  Use  (u(X),X)  as  the  new  (Uq,Xq)  and  go  back  to  Step  1. 

Note  that  very  often  the  computation  of  the  slope  u^  does  not  cause  much 
computational  overhead  because  we  usually  have  the  factorization  of  the 
Jacobian  computed  already  in  the  Newton  step.  Using  the  slope  in  a 
predictor-corrector  fashion  will  often  allow  us  to  take  a  much  bigger  step  in 
X  and  thus  reduce  the  overall  cost  of  determining  the  dependence  of  u  on  X. 

Unfortunately,  this  procedure  needs  some  modification  in  order  to  handle 
general  nonlinear  systems  because  of  the  possibility  of  existence  of  nonunique 
solutions.  The  nonuniqueness  usually  manifests  itself  in  the  form  of 
existence  of  'singular'  points  where  the  Jacobian  is  singular  (see  Figure 
3-1).  Points  such  as  point  A  in  Figure  3-1  are  called  limit  points  (or 
turning  points)  and  points  such  as  point  B  are  called  bifurcation  points. 
These  singular  points  are  further  characterized  by  the  conditions  that  £ 
Range (Gu)  at  a  limit  point  and  that  G^  e  Range (Gu)  at  a  bifurcation 
point  [12] • 
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The  difficulties  that  the  Natural  Continuation  Procedure  will  encounter 
are  three-fold.  First  of  all,  since  G^  is  singular  at  these  points,  Newton's 
method  will  in  general  be  only  linearly  convergent,  making  it  much  more  costly 
to  compute  the  solution.  Moreover,  when  tracing  a  solution  branch  near  a 
limit  point,  there  may  not  exist  a  solution  for  a  given  value  of  X  (see  Figure 
3-2)  and  hence  Equation  (3,5)  will  fail.  Lastly,  we  need  some  mechanism  for 
switching  branches  at  a  bifurcation  point. 

3.3  Arclength  Continuation 

In  the  pseudo-arclength  continuation  approach  [3,  12],  this  difficulty  is 
overcome  by  not  imposing  a  value  for  X  at  which  the  solution  u(X)  is  sought. 
Instead,  we  parametrize  the  solution  branches  using  an  arclength  parameter  s, 
and  specify  how  far  along  the  current  solution  branch  we  want  to  march. 

To  be  more  specific,  we  let  s  be  the  arclength  parameter,  and  treat  u($) 

•  « 

and  X(s)  as  functions  of  s.  Ve  can  compute  the  'slopes'  u(s^),  X(sq)  (where 
the  dots  denote  differentiation  with  respect  to  s)  of  a  known  solution  at  *s*q 
from  the  following  two  equations: 

e  • 

Gu  tt0  +  x0  GX  "  °*  (3'6) 

li;0ll2  +  U0|2  -  1  -  0.  (3.7) 

Equation  (3.6)  is  obtained  from  differentiating  G(u,X)  *  0  with  respect  to  s 
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and  Equation  (3.7)  imposes  tlie  arclength  condition.  We  can  theoretically 
'follow'  the  solution  curve  by  integrating  the  initial  value  problem  defined 

by  u(s)  and  X(s).  However,  this  process  may  become  unstable  -  the 

approximate  solution  to  the  initial  value  problem  may  deviate  more  and  more 
from  the  true  solution  curve.  One  way  to  stabilize  this  procedure  is  to  use 
u(s)  and  X(s)  in  a  predictor  step  to  provide  an  initial  guess  for  a 
Newton-type  method  to  bring  it  back  onto  the  true  solution  curve  (Figure  3-3). 

In  the  pseudo-arclength  continuation  procedure,  we  advance  from  s^  to  s 
along  the  solution  branch  by  requiring  the  new  solution  u($)  and  X(s)  to 
satisfy  the  following  equations: 

N(a (s).X(s))  -  -  a (sQ))  +  -  X(sQ>)  -  (s  -  sQ)  =  0 

and  (3.8) 

G(u( s) ,X( s) )  -  0.  (3.9) 

Equation  (3.8)  is  the  linearization  of  Equation  (3.7)  and  essentially  forces 
the  new  solution  to  lie  on  a  hyperplane  perpendicular  to  the  tangent  vector  of 
the  solution  curve  at  s^  and  at  a  distance  (s-Sq)  from  it.  Equation  (3.9) 
requires  u(s)  and  X(s)  to  lie  on  the  true  solution  curve  (Figure  3-3).  We  now 
solve  the  coupled  system  (3.8)  and  (3*9)  for  u(s)  and  X(s),  given  the  step 
size  (*~*q)  (efficient  strategies  for  choosing  the  step  size  are  discussed 
in  [21]).  We  can  use  Newton's  method,  in  which  case  we  have  to  solve  the 
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Figure  3-3:  Pseudo-Arclength  Continuation 


solution  curve  on  which 

G(u(a),A(j))  ■  0 


plane  J_  to  tangent  on  which 
N(u(o),X(o))  *  0 
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following  linear  system  at  each  iteration: 

T&uT  Tg  g  ~  ~suT  T  g  T 
I  I  I  u  x  I  II 

A  I  |*|  I  -  -I  I  (3.10) 

l6Xl  |NT  N,  8)J  I  N  I 

_  _  U  __  _  _  _ 

It  can  be  shown  that  at  limit  points,  where  G^  is  singular  and  i 
Range(G^),  the  linear  system  in  Equation  (3.10)  is  nonsingular  (see  [12])  and 
therefore  Newton's  method  for  the  coupled  system  (3.8)  and  (3.9)  is 
well-defined.  Hence  limit  points  present  no  problem  and  even  quadratic 
convergence  is  achievable. 

At  bifurcation  points,  where  G  is  singular  and  G1  e  Range (G  ),  things 

are  more  complicated.  In  the  simplest  case  of  only  one  branch  bifurcating 

from  the  main  branch  (simple  bifurcation),  an  additional  higher  order 

condition  (see  for  example  Crandall  and  Rabinowitz  [8])  involving  G  #  and 

G^  has  to  be  satisfied.  It  can  be  shown  [12]  that  this  condition,  together 

with  Equations  (3.6)  and  (3.7)  and  the  left  and  right  null  vectors  of  G^, 

.  • 

enable  two  solutions  for  (no#^o)  to  be  computed  at  a  simple  bifurcation  point, 

with  one  solution  corresponding  to  each  branch.  Using  the  appropriate  pair  of 
•  • 

U0,X0)  in  Equation  (3.8)  allows  branch  switching. 

In  order  to  solve  the  linear  system  in  Equation  (3.10)  by  direct  methods. 
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several  approaches  are  possible.  One  way  is  to  perform  Gaussian  elimination 

on  the  inflated  matrix  A,  with  some  form  of  pivoting  to  insure  stability.  But 
this  approach  completely  ignores  the  sparse  structure  which  is  usually  found 

in  Gu's  arising  from  nonlinear  elliptic  eigenvalue  problems.  In  order  to  take 
advantage  of  the  structure  in  G^,  K.ller  [12]  suggested  the  following 
block-elimination  procedure: 

Algpr it^  BE:  (block-elimination) 


Solve 

G 

u 

y 

-  Gk 

(3.11) 

and 

G 

u 

z 

-  -  G. 

(3.12) 

Set 

6k 

* 

1 

25 

W 

5E 

1 

N 

V 

1 

(3.13) 

and 

6u 

s 

z  -  6X  y. 

(3.14) 

Note  that  only  systems  with  the  coefficient  matrix  have  to  be  solved,  so 
structures  in  G^  can  be  exploited.  Moreover,  only  one  factorization  of  Gu  is 
needed.  It  has  been  shown  [23]  that  even  when  Gq  is  becoming  singular. 
Algorithm  BE  produces  iterates  that  converge  quadratically  at  limit  points. 

Continuation  methods  of  various  forms  and  levels  of  sophistication  have 
been  widely  used  in  the  engineering  literature  (see  the  survey  [3]).  For  a 
recent  survey  of  numerical  methods  for  bifurcation  problems,  see  for  example 
the  paper  by  Mittelmann  and  Weber  in  [19].  The  approach  taken  here  is  due  to 
Keller  [12],  and  has  recently  been  applied  to  other  problems  in  fluid 
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mechanics  (  [7],  [13]  ,  [15],  [22] ,  [23]  ).  A  related  approach  suggested 

by  Abbott  [1]  corresponds  (in  a  loose  way)  to  applying  Algorithm  BE  to  the 
matrix  A  with  the  last  column  permuted  into  the  first  n  columns  so  that  the 

corresponding  coefficient  matrix  in  Equations  (3.11)  and  (3.12)  becomes 

nonsingular  even  at  limit  points.  However,  as  has  already  been  pointed  out, 

any  structure  or  symmetry  in  is  lost  in  the  process,  and  hence  that 

approach  seems  unsuitable  for  large  elliptic  systems  in  two  or  three 

dimensions. 
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4.  Numerical  Implementation  for  the  Internal  Ware  Problem 

4.1  Computational  Domain  and  Boundary  Condition 

The  general  domain  of  interest  is  the  infinite  strip: 

Do  *  (  ~h2  <  y  <  H1  '  <  x  <  *  >  * 

In  the  numerical  implementation/  the  domain  is  truncated  to  a  rectangle  of 
large  but  finite  horizontal  extent.  If  the  solution  obtained  for  such  a  'box* 
is  to  represent  a  solitary  wave*  then  it  should  not  be  altered  if  the  domain 
is  enlarged.  Unless  otherwise  noted/  this  was  checked  numerically  for  all 
solutions  calculated.  The  truncated  domain  is  given  by: 

dl  =  (  ~h2  <  y  <  E1  '  *  L  <  x  <  L  } * 

and  the  boundary  condition  is  ; 

u  *  0  on  3D^. 

If  the  density  stratif ication  has  certain  symmetry/  then  the  domain  can  be 
further  reduced.  Specifically,  for  the  'tanhr  density  profile  (i.e.  F(y)  * 
tanh(y)  ) ,  F’(y)  is  symmetric  about  y  »  0,  and  it  can  be  shown  [24]  that 
antisymetric  solutions  satisfying  u(x,y,X)  -  -  u(x,-y,X)  exist.  Moreover, 
since  the  governing  equation  is  invariant  under  a  change  in  sign  in  the 
horizontal  coordinate  x,  solutions  symmetric  about  x  =  0  exist.  For  such  a 
wave,  the  domain  can  be  reduced  to  (Figure  4-1) 
D-(0<y<H,0<x<L)  , 
and  the  boundary  conditions  become 
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u  *  0  for  y  =  0  and  H  #  0  <  x  <  L 
and  for  x  *  L  ,  0  <  y  <  H  # 

dn/Bx  s  0  for  x*0,0<y<H. 

In  this  paper#  the  larger  domain  was  used  so  as  to  allow  for 

asymmetric  solutions  which  may  exist#  although  none  have  been  found  so  far. 

4.2  Computational  Mesh 

Since#  in  most  cases#  most  of  the  interesting  phenomenon  occurs  near  the 
origin#  it  would  be  a  waste  of  computational  effort  to  use  a  uniform  grid  on 
the  domain  D^.  The  approach  we  have  taken  is  to  transform  the  domain  into 
another  rectangular  domain  by  the  following  transformation  (Figure  4*1)  : 

5  =  tanh  (ox) 
tj  =*  tanh  (  0  y  )  . 

A  uniform  mesh  in  the  (g#n)  plane  would  put  more  mesh  points  near  the  origin 
of  the  domain  in  the  (x#y)  plane.  The  relative  spacing  of  the  mesh  points  in 
the  (x#y)  plane  can  be  varied  by  varying  the  parameters  a  and  0.  Note  that 
the  inverse  transformation  can  be  obtained  from  the  following  formula  : 

tanh”1  x  ■  (  log  (l+x)/(l-x)  )  /  2  if  lx I  <  1. 

Equation  (2.3)  now  becomes  : 

o2(i-?2)  a/a?  <i-?2)  a/a?  u 
+  p2(i-n2)  a/ati  (i-n2)  a/dn  a 


(4.1) 


rv 
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Figure  4-1:  Domain,  Boundary  Conditions  and  Mesh 
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+  X  u  sech2  (u  +  tanh-1  (n)  /  0)  =  0  . 


4.3  Discretization  and  Solution 

To  discretize  Equation  (4.1),  we  use  a  centered  second  order  finite 

difference  approximation.  Let  n  and  m  denote  the  number  of  intervals  in  the  £ 

and  r\  directions,  respectively,  and  let 

h  -  2  tanh(a  L)  /  n  ,  %  =  ^  +  ( i  -  1)  h  , 

k  =  2  tanh(p  H)  /  m  ,  +  (j  "  1)  k  . 

The  discretized  Equation  (4.1)  is  given  by  : 

(G(u,k)).  .  =  CE.  u .  -  .  +  CW.  n.  .  . 

i#J  i  x+l*j  i  i-l  ,j 


+  ui.j+i  +  csj  ui.j-i 

+  (  C5i +  Cl,j  >  tti.j 

2  -1 

+  A.  u.  .  sech  (u.  .  +  tanh  (q.)  /  0) 
* » J  i  >  J  j 


(4.25 


for  2  <=  i  <=  n  ,  2  <=  j  <*■  m  ,  where 

CE.  =  (1-52)  (l-^+1/2)  a2  /  h2 
CW.  =  (l-?2)  U-$i_1/2>  a2  / 

CN.  =  (1-t,2)  (l-t|j+1/2>  P2  /  *2  (4.3) 

cw.  *  (i-?2)  (i-52-1/2>  P2  7  k2 

C?i  "  -  (1-??+l/2  +  «-«i>  «2  '  h2 

C^j  *  -  (l-Hj+1/2  +  l-nJ.1/2>  <1-Wj)  P2  /  k2  • 
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The  derivatives  G^  and  needed  in  the  continuation  procedures  can  then 

be  easily  computed  from  Equation  (4.2).  The  discrete  operator  G  is  a 

u 

block-tridiagonal  matrix  of  sixe  (n-l)*(m-l)  and  has  the  familiar  nonzero 
structure  generated  by  five-point  operators.  Since  this  matrix  is 
nonseparable ,  we  cannot  use  fast  Poisson  solvers.  Moreover,  the  matrix  is 
indefinite  and  hence  iterative  methods  like  Successive-Over-Relaxation  cannot 
be  used  directly.  We  have  chosen  to  use  a  band  solver  called  LEQT1B  from 
IMSL  [11].  This  choice  was  satisfory  for  our  problem  in  terms  of  speed  and 
storage.  Typically,  L  =  20  was  found  to  be  large  enough  to  contain  the 
solitary  waves.  Various  values  of  the  half-depth  H  were  used,  ranging  from  H 
=4  to  H  =  40.  Results  from  the  H  -  40  case  was  meant  to  be  compared  to 
experimental  results  obtained  by  Davis  and  Acrivos  [9].  Typically,  n  =  20  and 
m  =  32  was  used.  These  were  found  to  provide  adequate  resolution.  All 
computations  were  performed  on  a  GDC  Cyber  176  computer,  and  it  takes  slightly 
less  than  one  CPU  second  to  obtain  a  solution  u  for  a  given  value  of  X. 
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5.  Numerical  Results 

In  this  section,  we  present  a  sample  of  the  numerical  results  that  we 
have  obtained  for  the  internal  wave  problem  with  the  'tanh'  density  profile, 
i.e.  F'(y)  =  sech2(y) .  All  waves  presented  here  are  'mode-2'  waves 
(anti-symmetric  about  y  =  0)  and  are  symmetric  about  x  *  0.  Other  waves 
corresponding  to  different  density  profiles  and  different  modes  have  been 
computed,  but  won't  be  presented  here  because  they  don't  correspond  to  the 
experimental  results  of  Davis  and  Acrivos  [9]  and  probably  are  not  as  easily 
realizable  as  the  mode-2  waves.  More  results  and  discussions  on  the  fluid 
dynamical  significance  of  these  computed  waves  can  be  found  in  [24] • 

The  dependence  of  u  on  X  as  computed  is  shown  in  Figure  5-1  for  H  »  4,  10 
and  40.  Each  value  along  these  curves  is  a  nonlinear  solution  corresponding 
to  Equation  (4.2),  and  many  such  solutions  (15  to  20)  are  used  to  construct 
each  curve. 

For  the  shallow  water  case  of  H  =  4,  we  have  found  nonunique  solutions 

and  a  limit  point  at  X  *  0.979114  with  (-u)  =  2.647.  Some  contour  plots  of 

max 

the  total  stream  function  u+y  for  various  values  of  X  are  presented  in  Figures 
5-2,  5-3  and  5-4.  Contour  plots  for  the  case  H  *  40  are  presented  in  Figures 
5-5,  5-6  and  5-7. 


Horizontal  cross-sections  cut  along  a  line  passing  through  the  y-location 
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of  the  wave  maximum  are  displayed  in  Figure  5-8.  The  vertical  cross-sections 
through  the  same  maximum  are  displayed  in  Figure  5-9.  The  uppermost  line  in 
Figure  5-8  is  a  portion  of  the  extremely  large  amplitude#  bore-like  wave  found 
at  X  *  1.1615,  past  the  limit  point.  This  wave  probably  does  not  correspond 
to  a  physical  solitary  wave  because  the  artificial  boundary  seems  to  have  a 
significant  influence  on  it.  However,  it  does  show  that  the  mathematical 
problem  given  by  Equation  (4.2)  has  nonunique  solutions.  Ve  should  note  here 
that  we  have  not  addressed  the  question  of  stability  of  these  steady  state 
solutions. 

In  Figure  5-10,  the  experimental  results  of  Davis  and  Acrivos  [9]  are 
compared  with  the  computed  results.  The  density  profile  in  their  experiment 
can  be  well  approximated  by  F(y)  =  tanh(y)  and  the  half-depth  is  equivalent  to 
our  case  of  H  =  40.  The  relative  density  change  across  the  stratified  layer 
in  all  the  experiments  are  relatively  small  (i.e.  a  is  small  in  Equation 
(2*1)),  and  hence  the  Boussinesq  approximation  made  to  arrive  at  Equation 
(2.3)  is  justified.  Also  displayed  in  Figure  5-10  is  the  analytical  result 
from  the  weakly  nonlinear  theory.  As  can  be  seen  from  the  figure,  excellent 
agreement  is  obtained  with  the  fully  nonlinear  results  of  our  study,  even  in 
the  large  amplitude  and  nonlinear  regime,  where  the  weakly  nonlinear  theory 
starts  to  deviate  from  observed  data. 


Figurs  5-5:  Contours  of  n  for 


Contour 


Figure  5-10:  Cooper i ton  with  Date  of  Dawis  end  Acrivos 


6 .  Cone las ion 


In  this  paper,  we  described  efficient  continuation  techniques  for  solving 
general  nonlinear  eigenvalue  problems.  We  applied  these  techniques  to  the 
solution  of  Long's  Equation  which  governs  the  propagation  of  internal  solitary 
waves  in  a  stratified  medium.  Nonunique  solutions  have  been  found  by  the 
continuation  procedures.  Computed  results  show  excellent  agreement  with  the 
experimental  data  obtained  by  Davis  and  Acrivos  [9],  even  in  the  nonlinear 
regime  where  weakly  nonlinear  theories  fail  to  give  good  predictions. 
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